
# Estimating the persistence of party cue effects in a panel survey experiment

# SENSITIVITY POWER ANALYSIS FILE ####

library(tidyverse)
library(lme4)
library(broom.mixed)

a <- 0   # population average intercept
b <- 0.2 # population average treatment effect

sigma_a_sub <- 0.4 # std dev in intercepts across subjects
sigma_b_sub <- 0.2 # std dev in treatment effects across subjects

sigma_a_issue <- 0.3 # std dev in intercepts across issues
sigma_b_issue <- 0.1 # std dev in treatment effects across issues

sigma <- 1 # std dev in outcome

# Subject and issue Ns
n_subjects <- 773
n_issues   <- 5

# Write function to simulate data and fit model
fun_run_sims <- function(run) {

  # Draw intercepts and treatment effects for each subject and issue
  vary_fx_subjects <- 
    data.frame(a_sub = rnorm(n = n_subjects, mean = a, sd = sigma_a_sub),
               b_sub = rnorm(n = n_subjects, mean = b, sd = sigma_b_sub),
               subject = 1:n_subjects)
  
  vary_fx_issues <- 
    data.frame(a_issue = rnorm(n = n_issues, mean = a, sd = sigma_a_issue),
               b_issue = rnorm(n = n_issues, mean = b, sd = sigma_b_issue),
               issue = 1:n_issues)
  
  # Bind into data frame, generate treatment assignments and draw opinion data
  d <- 
    expand_grid(vary_fx_subjects, vary_fx_issues) %>% 
    mutate(cue = rbinom(nrow(.), 1, 0.5)) %>% # randomize treatment assignment
    mutate(mu = a + a_sub + a_issue + b*cue + (b_sub - b)*cue + (b_issue - b)*cue) %>% # compute mu
    mutate(opinion = rnorm(n = n(), mean = mu, sd = sigma)) # draw observed opinions
  
  # Fit model
  fit <- lmer(formula = opinion ~ 1 + cue + (1 + cue || issue) + (1 + cue || subject), data = d)
  
  # Tidy and index
  out <- broom.mixed::tidy(fit) %>% 
    mutate(singular = lme4::isSingular(fit),
           run = run)
  
  out
  
}

# How many simulations?
n_sim <- 2000

set.seed(42)

# Run simulations
sim_out <- map(as.list(c(1:n_sim)), ~fun_run_sims(run = .x))

# Analyze results
sim_df <- bind_rows(sim_out)
saveRDS(sim_df, "power_sims_df.rds")
# sim_df <- readRDS("power_sims_df.rds")

critval <- qnorm(.025, lower.tail = F)

sim_df <-
  sim_df %>% 
  mutate(p_less_than_05 = case_when(statistic > critval ~ 1,
                                    statistic < critval ~ 0)) %>% 
  mutate(term_c = case_when(effect == "ran_pars" ~ str_c(group, term, sep = "_"),
                            effect == "fixed" ~ term)) %>% 
  mutate(term_c = str_remove(term_c, ".1"))

# How many lme4 fits were singular?
sim_df %>% 
  group_by(run) %>% 
  sample_n(1) %>% 
  ungroup() %>% 
  count(singular)

sim_df %>%
  filter(singular == FALSE) %>% 
  group_by(term_c) %>% 
  summarise(mean_value = mean(estimate)) # simulation recovers the true parameter values

sim_df %>% 
  filter(singular == FALSE) %>% 
  filter(term == "cue") %>% 
  summarise(prop_p = mean(p_less_than_05)) # power estimate on the fixed ATE

